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Introduction 

This  is  the  second  installment  of  a  report-pair  concern- 
ing implementation  of  tensor  product  factoring  of  coeffi- 
cient matrices  in  applications  of  the  finite  element  method 
to  numerical  weather  prediction.  It  was  noted  in  Part  1 
(Ref.  1)  that  these  techniques  were  introduced  in  numerical 
weather  prediction  by  Staniforth  and  Mitchell  (Ref.  2). 
Discussed  in  Part  1  are  applications  in  which  the  "mass" 
matrix  for  a  grid  such  as  that  shown  in  Fig.  1  is  factored 
as  the  tensor  product  of  two  matrices. 
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Fig.  1.   Node  numbering  and  spacing. 

One  of  these  matrices  (MA)  depends  solely  upon  the  nodal 
spacing  in  the  east-west  direction  (a.)  and  the  other  (MB) 
depends  only  on  the  north- south  spacing  (b-)«  We  began  with 
the  set  of  simultaneous  linear  equations 

M  w  =  v,  <1> 

where  M  (the  "mass"  matrix)  is  symmetric,  ne  x  ne ,  and  w  and 
v  are  column  vectors  of  height  ne .  M  and  v  are  input  quan- 
tities and  w  is  sought.  The  tensor  product  representation 
of  M  is 

M  =  MB  *  MA,  <2> 

where  MB  and  MA  are  tridiagonal,  symmetric  matrices,  e  x  e 
and  n  x  n,  respectively.  (The  tensor  product  and  matrices 
MA  and  MB  are  defined  in  Appendix  A. )  This  representation 
allowed  <1>  to  be  rewritten  as 

MA  W  MB  =  V,  <3> 

where  W  is   n  x  e  and  the  successive   columns  are  subvectors 

3 


of  w  corresponding  to  the  rows  of  Fig.  1.  V  is  also  n  x  e 
and  similarly  derived  from  v.  Boundary  conditions  consid- 
ered were  a  cyclic  condition  in  the  east-west  direction  and 
either  homogeneous  Neumann  conditions  (normal  derivative 
zero)  or  nonhomogeneous  Dirichlet  conditions  (specified 
nonzero  values)  on  the  northern  and  southern  edges. 

The  present  report  discards  the  cyclic  east-west  boundary 
condition  and  deals  with  two  cases: 

(1)  Solutions  of  <3>  with  nonhomogeneous  Dirichlet  con- 
ditions on  all  four  edges; 

(2)  Solution  of  Poisson's  equation  for  the  same  region 
with  nonhomogeneous  Dirichlet  conditions  on  all  four 
edges . 

Mass  Matrix  ^_   Dirichlet  Boundary  Conditions 

Effects  of  the  Dirichlet  boundary  conditions  on  the  solu- 
tion process  are  most  readily  understood  by  considering  the 
following  partitioned  form  of  <1>: 

Dfeiftsiiyg-M      '    <4> 

In  <4>  the  w  vector  has  been  rearranged  so  that  all  of  the 
boundary  values  are  in  the  subvector  wu  and  the  interior 
("center")  values  are  in  wc .  A  similar  reordering  has  been 
applied  to  v  and  M.  If  the  boundary  values  of  w  are  pre- 
scribed, then  w^  is  known  and  only  wc  remains  to  be  found. 
Expanding  the  lower  partition  of  <4>  and  placing  the  known 
terms  on  the  right  gives 

M22WC  =  vc  ~  M2iWb»  <5> 

or,  letting  v  '  =  v   -  M21W,  ,  we  have 
°   c     c       b 

M22wc  =  vc' .  <5f  > 

We  consider  now  how  the  strategy  just  described  can  be 
applied  when  the  tensor  product  resolution  of  M  has  been 
used  to  convert  <1>  into  <3>.  In  the  matrix  W  the  pre- 
scribed boundary  values  occupy  the  first  and  last  columns 
and  the  top  and  bottom  rows.  Denote  this  border  matrix, 
including  an  (n-2)  x  (e-2)  null  matrix  inside,  by  WB . 


Calculate 

VB  =  MA  WB  MB  <6> 

and  now  form 

V  =  V  -  VB.  <7> 

Now  define  a  set  of  submatrices  MAI,  MB1,  Wl,  and  VI 
obtained  from  MA,  MB,  W,  and  V',  respectively,  by  removing 
the  first  and  last  columns  and  the  top  and  bottom  rows.  The 
reduced  problem  becomes 

MAI  Wl  MB1  =  VI  <8> 

As  described  in  Ref.  1,  <8>  may  be  solved  by  standard  Gaus- 
sian elimination  procedures.  A  computer  program  (GAUSS4) 
which  carries  out  these  calculations  is  listed  in  Appendix 
B.  The  subroutines  of  GAUSS4  are  designed  for  substitution 
in  the  program  devised  by  Hinsman  (Ref.  5). 

Poisson'  s  Equation  -_  Dirichlet  Boundary  Conditions 

As  noted  above,  Stanif orth (and  Mitchell  (Ref.  2)  appear 
to  have  been  first  in  applying  the  tensor  product  resolution 
to  Poisson' s  equation  in  a  numerical  weather  prediction 
problem  using  the  finite  element  method.  Additional  detail 
is  given  in  earlier  papers  by  Dorr  (Ref.  3)  and  by  Lynch, 
Rice,  and  Thomas  (Ref.  4). 

Finite  element  discretization  of  Poisson' s  equation  for 
the  region  of  Fig.  1  results  in  a  set  of  simultaneous  linear 
equations  which  may  be  written  in  matrix  form  as 

K  w  =  v,  <9> 

where  vectors  v  and  w  are,  respectively,  given  and  unknown. 
As  for  <1>,  each  has  length  ne  and  the  coefficient  matrix  K 
is  ne  x  ne ,  symmetric,  sparse,  and  block  tridiagonal.  K  is 
called  the  "stiffness"  matrix  in  finite  element  parlance. 

It  is  easily  shown  that  K  is  expressible  as  the  sum  of 
two  tensor  products  as  follows: 

K  =  MB  *  SA  +  SB  *  MA.  <10> 

The  new  matrices  SA  and  SB  are  symmetric,  tridiagonal  and 
depend  only  on  the  a.  and  b.,  respectively.  Explicit  formu- 
las for  SA  and  SB  are  given  in  Appendix  A. 
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Using  the  definition  of  the  tensor  product  and  again  con- 
verting the  vectors  w  and  v  into  the  n  x  e  rectangular 
matrices  W  and  V,  <9>  may  be  written  as 

SA  W  MB  +  MA  W  SB  =  V.  <11> 

Before  solving  <11>  we  must  first  take  account  of  the  Diri- 
chlet  boundary  conditions  on  the  four  edges  of  the  region. 
As  in  solving  <3>,  the  given  boundary  values  are  in  the 
first  and  last  columns  and  top  and  bottom  rows  of  W.  As 
before,  we  let  WB  be  an  n  x  e  matrix  containing  the  given 
boundary  values,  together  with  zeros  at  locations  corre- 
sponding to  interior  nodes.   Calculate 

VB  =  SA  WB  MB  +  MA  WB  SB,        <12> 
and  then  form 

V  =  V  -  VB.  <7> 

The  remaining  step  again  parallels  that  used  when  applying 
the  Dirichlet  boundary  conditions  to  <3>.  Specifically,  we 
introduce  submatrices  MAI,  MB1,  SA1 ,  SB1,  Wl,  and  VI 
obtained  from  MA,  MB,  SA,  SB,  W,  and  V',  respectively,  by 
removing  the  first  and  last  columns  and  the  top  and  bottom 
rows.   The  reduced  problem  becomes 

SA1  Wl  MB1  +  MAI  Wl  SB1  =  VI.    <13> 

To  solve  <13>  we  first  need  the  complete  solution  of  the 
eigenproblem 

SB1  p.  =  Xi  MB1  Pl,  <14> 

where  p-  is  the  ith  eigenvector  and  \-  is  the  corresponding 
eigenvalue.   We  write  the  complete  solution  in  the  form 

SB1  P  =  MB1  P  A,  <14'> 

where  P  is  the  (e-2)  x  (e-2)  modal  matrix  whose  columns  are 
the  p.  and  A  is  the  (diagonal)  spectral  matrix  whose  ele- 
ments are  the  X,*  .  We  specify  that  the  modal  matrix  is  nor- 
malized so  that 

PT  MB1  P  =  I,  <15> 

where  I  is  the  identity  matrix  of  order  e-2  and  PT  is  the 
transpose  of  P.  If  both  sides  of  <13>  are  postmultiplied  by 
P  and  <14'>  is  used  to  replace  SB1  P,  <13>  becomes 

SA1  Wl  MB1  P  +  MAI  Wl  MB1  PA  =  VI  P.      <16> 
Let  X  =  Wl  MB1  P  and  U  =  VI  P,  then  <16>  is  equivalent  to 
(SA1  +  X-  MAI)  x.  =  u.,     1=1,  e-2,     <17> 
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where  x-  and  u-  are,  respectively,  the  ith  columns  of  X  and 
U.  Since  the  coefficient  matrix  in  <17>  is  tridiagonal ,  the 
Gsiussian  elimination  process,  i.e.,  factoring,  forward 
reduction,  and  back- substitution,  is  computationally  econom- 
ical. The  final  step  consists  of  a  matrix  multiplication  to 
obtain 

Wl  =  X  PT.  <18> 

Since  Wl  contains  the  w  values  at  all  interior  nodes  and  the 
boundary  values  were  known  in  advance,  the  solution  is  com- 
plete. A  FORTRAN  program  (GAUSS5)  which  implements  the  ten- 
sor product  solution  for  Poisson's  equation  is  given  in 
Appendix  C. 

Operation  Counts  and  Storage  Requirements  2.  Poisson' s  Equa- 
tion 

In  Ref.  1  comparisons  of  floating  point  operation  counts 
and  storage  requirements  were  made  for  solutions  of  <1>. 
Substitution  of  the  boundary  conditions  considered  here  in 
place  of  those  considered  in  Ref.  1  has  a  negligible  effect 
on  both  operation  counts  and  storage  requirements.  Accord- 
ingly, no  further  comparison  is  given  here  for  solutions  of 
<1>. 

Solution  of  Poisson's  equation  (<9>)  using  the  tensor 
product  resolution  <10>  of  K  is  more  costly  in  terms  of  com- 
putation and  storage  than  the  previously  studied  applica- 
tions to  <1>.  In  Table  1  the  number  of  floating  point  oper- 
ations and  the  required  number  of  coefficient  matrix  storage 
locations  are  compared  for  three  different  solution  methods. 
These  are  SOR  (successive  over-relaxation),  SKY  (skyline 
storage  and  Gauss  elimination) ,  and  TENSOR  (the  scheme 
described  above).  A  floating  point  operation  is  defined  to 
be  one  multiplication  (or  division)  plus  one  addition  (or 
subtraction) .  The  exact  operation  counts  would  be  polynomi- 
als in  n  and  e.  Only  the  highest  degree  terms  are  given  in 
the  table.  Since  it  is  not  possible  to  predict  the  number 
of  iterations  per  solution  using  SOR,  the  operation  count 
given  for   that  algorithm   is  for   a  single   iteration.    In 
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Table  1  a  storage  location  corresponds  to  8  bytes, 
comparison  it  is  assumed  that  each  floating  point  number 
requires  8  bytes  of  storage  and  an  integer  requires  4  bytes. 
The  storage  requirement  given  for  SOR  is  based  on  the  com- 
pact  storage    scheme   described   by   Franke    and   Salinas 

(Ref.  6). 

► 

TABLE  1.   Operation  Counts  and  Storage  Requirements. 


ALGORITHM 

NUMBER  OF  OPERATIONS 

NUMBER  OF  STORAGE  LOCATIONS 

PER  SOLUTION 

FOR  COEFFICIENT  MATRICES 

SOR 

10  en     (1) 

13  en 

SKY 

2  en2 

en2 

TENSOR 

2  en2 

e2 

Note:  1.  Number  of  operations  per  iteration. 

It  is  perhaps  surprising  to  note  that  the  number  of  oper- 
ations for  TENSOR  is  no  fewer  than  for  SKY.  Turning  atten- 
tion to  storage  requirements  reveals  that  for  a  large  prob- 
lem (e  =  n  =  100,  say)  the  SKY  storage  requirement  for  the 
stiffness  matrix  is  8  megabytes,  compared  with  1  megabyte 
for  SOR  and  80  kilobytes  for  TENSOR.  It  is  this  comparison 
which  is  the  compelling  reason  for  preferring  TENSOR.  It  is 
acknowledged  that  there  is  overhead  associated  with  the  one- 
time solution  of  the  eigenvalue  problem  <14>,  but  the  tri- 
diagonal  form  of  matrices  SB1  and  MB1  makes  the  amount  of 
computation  comparable  with  that  required  for  a  single  solu- 
tion of  Poisson's  equation.  Since  two  solutions  of  Pois- 
son's  equation  are  required  at  each  time  step,  the  overhead 
is  clearly  negligible. 

It  is  not  feasible  to  make  a  definitive  comparison 
between  the  number  of  operations  required  for  SOR  and  those 
required  for  the  other  two  algorithms.  If  the  number  of 
iterations  is  less  than  0.2  e,  then  SOR  will  be  more  econom- 
ical and  the  storage  tradeoff  would  need  to  be  weighed. 


Conclusions 

It  has  been  demonstrated  that  Dirichlet  boundary  condi- 
tions on  all  edges  of  the  region  are  easily  incorporated  in 
solution  processes  which  use  tensor  product  resolution  of 
the  coefficient  matrix.  For  very  large  problems  the  tensor 
product  algorithm  uses  much  less  core  storage  than  alterna- 
tive choices.  The  computational  expense  of  a  solution  to 
Poisson's  equation  is  substantially  the  same  for  Gaussian 
elimination  and  for  the  tensor  product  scheme.  It  is 
expected  that  successive  over-relaxation  is  almost  always 
more  expensive. 
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APPENDIX  A  -  TENSOR  PRODUCT  AND  MATRIX  DEFINITIONS 


Tensor  Product 


The  tensor  product   of  matrices  C   and  D 


may  be  represented  in  block  partition  form  as 


C*D  = 


Cn  D   c12D   Ci3D" 

c2i  D   c22  D   c23  D 

.C31O   c32D   C33D 


where  the  c-jj  are  the  elements  of  C.  Note  that,  if  C  and  D 
have  dimensions  r  x  s  and  t  x  u,  respectively,  the  tensor 
product  has  dimensions  rt  x  su . 

Definitions  for  matrices  MA  and  SA  are  given  below.  The 
corresponding  expressions  for  MB  and  SB  may  be  obtained  by 
substituting  "b"  for  "a"  throughout  and  replacing  n  by  e . 
(Symbols  a-j  and  b -j  are  defined  in  Fig.  1.) 
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APPENDIX  B 

PROGRAM  to  SOLVE  M  vr  =  V  with  DIRICHLET  EOUNDARY  CONDITIONS 

Listing:   GAUSS4  FORTRAN 

C   MAIN  PROGRAM   .MASS  MATRIX  USING  TENSOR  PRODUCT  RESOLUTION 

C 

C   THIS  PROGRAM  IS  DESIGNED  TO  TEST  THE  SCHEME  (TENSOR) 

C   WHICH  RESOLVES  THE  MASS  MATRIX  INTO  A  TENSOR  PRODUCT  IN 

C   ORDER  TO  SOLVE  THE  SYSTEM  OF  EQUATIONS  M  w  =  v . IN 

C   THIS  PROGRAM  THERE  ARE  DIRICHLET  BOUNDARY  CONDITIONS  ON 

C   ALL  4  EDGES  OF  THE  REGION.   THE  PRESCRIBED  BOUNDARY 

C   VALUES  ARE  GIVEN  IN  THE  CORRESPONDING  LOCATIONS  IN  V. 

C   THE  SUBROUTINES  MAY  BE  INSERTED  IN  THE  PROGRAM  DEVISED 

C   BY  HINSMAN. 

C 

IMPLICIT  REAL*8(A-H.O-Z) 

COMMON/ CM 1A/ NLAT, NLONG 

COMMON/ CM8/A(Z1KB7z1] 

COMMON  AG(ZB),BG(ZC) ,GAD(ZK) ,GBD(ZL) ,MA(ZM) ,MB(ZN) 

DIMENSION  V(ZP) 

READ(5,  *) NLONG,  NLAT 

LATX=NLAT+1 

WRITE (6, 1000) 
1000   FORMAT (//,'    MASS  MATRIX  -  TENSOR  PRODUCT  RESOLUTION 


T  I 


s& 


TE ( 6 , 1 0  0 1 ) NLONG , NLAT 


READ ( 5  *) A. B 

WRITE(6,500)A 
WRITEI6.503)B 

503  FORMAT?/,'    B:   ',(24F3.0)) 

500  FORMATf/,'    A:   ,,[24F3.0)) 

1001  FORMAT (*     NLONG  =  ,13,'      NLAT  =  ',13  ,/J 
C   CONSTRUCT  FACTORS,  GAD  AND  GBD,  OF  MASS  MATRIX 

CALL  AMTRX3 
WRITE (6. 501) AG 

501  FORMATf )'         AG:   ' ,(12F4.1)) 
WRITE(6;504)BG 

504  FORMAT ()     '         BG :   ',(12F4.1)) 
WRITE?6,l002)GAD 

1002  FORMATS,  '   GAD  '  ,  /  ,  ( 3X ,  6F7  .  3  )  ) 
WRITE (6, 1004 )GBD 

1004   FORMAT?//.'   GBD' ,/ ,(3X,6F7.3)) 


WRITE(6,1006)MB 
K=?NLAT-1) -NLONG 
L= NLAT -NLONG 
READ(5  *)V 
WRITE? 6, 5 10 )V 
C   CORRECT  RIGHT-HAND  SIDE  FOR  DIRICHLET  CONDITION 
LONGM=NLONG-l 
DO  2  J=2.LONGM 

V ( J+NLONG ) =Y ( J+NLONG ) - (GAD ( 2 * J- 1 ) *V ( J- 1 ) +GAD ( 2 -J- 2 ) 
l*V(Jl  +  GAD(2*M)*V(J+lh*GBD(3) 

2  V(K+J)=V(K+J)-(GAD(2"J-1)"V(L+J-1)+GAD(2"J-2)"V(L+J) 
l+GADi2*J+l)*V(L+J+I))*GBD(2*LATX-l) 

CU=GBD(3) 

CL=GBD(2"LATX-1) 

GBD(3)=0. 

GBDf2--LATX-l)=0. 

DO  3  J =2, NLAT 

Y((J-l)"NLONG+2)=V( (J-l)-NLONG+2)- (GBD(2"J-1)-V( (J-2) 
1 -NLONG  + l) + GBD ( 2 -J - 2 1*V ( ( J- 1 ) "NLONG  + 1 ) + GBD ( 2 -J  + 1 ) 
2*vTj*NLONG+l)J*GAD(3) 

3  V(J"NLONG-l)=V(J"NLONG-l)- (GBD(2"J-1)"V( (J-l)-NLONG) 
1  +  GBD ( 2  - J- 2 ) -V ( J -NLONG ) +GBD ( 2  - J+ 1 ) *V ( ( J+ I ) -NLONG ) ) 
2*GADf2*NLONG-l) 

GBD(3)=CU 
GBD_f2"LATX-l>CL 
WRITE (6, 5 10 )V 
C   PERFORM  LDLT  FACTORING  OF  GAD  AND  GBD 

4  CALL  FACT 1( GAD, NLONG) 
CALL  FACTlfGBD.LATX) 
WRITE(6,10O2)GAD 
WRITE (6, 1004 )GBD 
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C   PERFORM  FORWARD  REDUCTION  AND  BACK- SUBSTITUTION  USING 
C   FACTORS  OF  GAD 

CALL  BACKA1(GAD,V) 


WRITE(6,510)V 
LFORM  FORWARD 


C   PERFORM  FORWARD  REDUCTION  AND  BACK- SUBSTITUTION  USING 
C   FACTORS  OF  GBD 

CALL  BACKBl(GBD.V) 
6      WRITeT6,510]V 

510    FORMATf/,'  V:  ' ,5F8 . 2 . / , (4X,5F8 .2) ) 
1003   FORMATf/,'    MA: *  , 2X ,3613 ) 
1006   FORMAT (/,'    MB : ' , 2X, 3613 ) 

STOP 

END 

SUBROUTINE  FACT1(A,NN) 
C 

C  SUBROUTINE  FACT1  PERFORMS  L*D*LT  FACTORING  ON  A  SUBMATRIX 

C  OF  A  SYMMETRIC  TRIDIAGONAL  MATRIX  STORED  IN  SKYLINE  FORM. 

C  THE  SUBMATRIX  IS  FORMED  BY  OMITTING  THE  FIRST  AND  LAST 

C  COLUMNS  AND  ROWS  OF  THE  INPUT  MATRIX. 

C  . -  -  INPUT  VARIABLES 

C  .    A(NWK)  =  INPUT  MATRIX  STORED  IN  COMPACTED  FORM 

C  .    NN      =  NUMBER  OF  COLUMNS  (OR  ROWS)  IN  INPUT  MATRIX 

C  .    NWK     =  NUMBER  OF  ELEMENTS  BELOW  SKYLINE  (2*NN  -  1) 

C  .   -  -  OUTPUT  -  - 

C  .    A (NWK)  =  D  AND  L  -  FACTORS  OF  INPUT  SUBMATRIX 

C  . 

C  . 


IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  A(l) 
C 
C      PERFORM  L*D*LT  FACTORIZATION  OF  STIFFNESS  MATRIX 


C 


LONGMM=NN-2 

A(3)=0. 

DO  50  J=2,LONGMM 


120 

STOP 
50     A?2*J+1)=TEMP 
2000   FORMATf//, '  STOP    MATRIX  NOT  POSITIVE  DEFINITE'  //, 

1'  NONPOSITIVE  PIVOT  FOR  EQUATION   ',14,//,'  PIVOT  =  ', 
2E20.12) 
RETURN 
END 

,,,„,„,SyBRpyTJNE  ,JA£l^nA,  V^    ^  ^   ,  ^ 

C   THIS  SUBROUTINE  PERFORMS  THE  FORWARD  REDUCTION  AND  BACK- 

C   SUBSTITUTION  USING  THE  FACTORS  OF  GAD 

C 

IMPLICIT  REAL-8(A-H.O-Z) 

COMMON/ CM 1A/NLAT .NLONG 

DIMENSION  A(1),V(1) 
C 

C      DEFINE  LIMITS  FOR  DO-LOOPS 
C 

NTM=NLAT-1 

LONGM=NLONG-l 

LONGMM=NLONG-2 
C 

C      REDUCE  RIGHT-HAND- SIDE  LOAD  VECTOR 
C 

DO  100  K=1,NTM 

DO  20  J=3,LONGM 
20     V (K*NLONG+ J ) =V (K-NLONG+ J ) -V (K-NLONG+ J- 1 ) *A( 2 * J- 1 ) 

C      DIVIDE  BY  DIAGONAL  ELEMENTS 
C 

DO  40  J=l,LONGMM 
40     V(K*NLONG+J+l)=V(K»NLONG+J+l)/A(2*J) 

C      BACK- SUBSTITUTE 
C 
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DO  60  J=3,LONGM 

L=(K+ll"NLONG-J+l 
M= 2 "(NLONG -J) +  3 

60   v{l)=v7l)-v(l+i)*a(m) 

100   continue 


PFTURN 


END 
C 

C   THIS  SUBROUTINE  PERFORMS  THE  FORWARD  REDUCTION  AND  BACK- 

C   SUBSTITUTION  USING  THE  FACTORS  OF  GBD . 

C 

IMPLICIT  REAL*8(A-H,0-Z) 

COMMON/CM 1A/ NLAT. NLONG 

DIMENSION  A(l) ,V(1) 
C 

C   DEFINE  NEEDED  INDEX  VARIABLES 
C 

LATX=NLAT+1 

LONGM=NLONG-l 
C 

C      REDUCE  RIGHT-HAND- SIDE  LOAD  VECTOR 
C 

DO  100  K=2,LONGM 

DO  20  J=3,NLAT 
20     V(K+  (J-l)-'NLONG)=V(K+  ( J-  1 )  -NLONG  )  -V(K+  (J-2)*NL0NG) 

1*A(2*J-1) 
C 

C   DIVIDE  BY  DIAGONAL  ELEMENTS 
C 

DO  40  J=2,NLAT 

C      BACK- SUBSTITUTE 
C 

DO  60  J=3,NLAT 
60     V(K+ (LATX-j7*NLONG)=V(K+ (LATX- J)-NLONG) 

1-V(K+(1+LATX-J)*NL0NG)*A(2*(LATX'-J)  +  3) 
100    CONTINUE 

RETURN 

C   THIS  SUBROUTINE  FORMS  THE  MASS  MATRIX  IN  THE  FORM  OF  A 

C   TENSOR  PRODUCT  OF  THE  GBD  MATRIX  AND  THE  GAD  MATRIX. 

C   THE  FIRST  OF  THESE  IS  NLAT  +  1  BY  NLAT  +  1,  SYMMETRIC, 

C   AND  TRIDIAGONAL.   THE  SECOND  IS  NLONG  BY  NLONG,  SYMMET- 

C   RIC,  AND  TRIDIAGONAL.   NOTE  THAT  THERE  ISNO  CYCLIC 

C   BOUNDARY  CONDITION  IN  THE  EAST-WEST  DIRECTION.   BOTH  GBD 

C   AND  GAD  ARE  STORED  IN  SKYLINE  VECTOR  FORM  ?UPPER  TRIANGLE 

C   WITH  SPACE  FOR  FILL-IN).  INTEGER  ADDRESS  VECTORS  MB  AND 

C   MA  ARE  ALSO  GENERATED. 

C 

IMPLICIT  REAL*8(A-H.O-Z) 

COMMON/ CM 1A/ NLAT , NLONG 

COMMON/ CM8/A(Z1) >B(Z1) 

COMMON  AG(ZB) .BGIZC) .GAD(ZK) , GBD(ZL) ,MA( ZM) ,MB(ZN) 
C      DIMENSION  BG ( NLAT ), AG (NLONG ),GBD(2"NLAT-1), 
C     lGAD(3*NLONG-3) ,MA(nLONG+1) ,MB(NLAT+2) 

LATX=NLAT+1 

LONGM=NLONG-l 
C    FIND  BG  =  [ELEMENT  HEIGHT) /6. 

LONGM=NLONG-l 

DO  2  J=1,NLAT 
2      BG(J)=B(l+LONGM*(J-l))/3. 
C    GENERATE  GBD 


GBD(1)=2.*BG(1) 
DO  4  J=2,NLAT 
K=2*(J-l) 
GBD(K)=2.*(BG(J-1)+BG(J)) 


4      GBD(K+1)=BG(J-1) 

GBD ( 2 -NLAT )  =  2 . -BG (NLAT ) 
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GBD ( 2*NLAT+  L )  =BG (NLAT ) 
C    FIND  AG  =  (ELEMENT  WIDTH)/ 6 

DO  10  J= I, LONGM 
10     AG(JI=A(JJ/3. 
C    GENERATE  GAD 

GAD(1)=2.*AG(1) 

DO  12  J= 2, LONGM 

K=2*(J-1) 

GAD(K)=2. *(AG(J-1)+AG(J) ) 
12     GAD  K+1)=AG(J-1) 

GAD  ( 2  *  LONGM )  =  2*AG  ( LONGM ) 
GAD( 2*LONGM+ 1} =AG f LONGM) 
C    GENERATE  DIRECTORY  VECTORS 
MB ( 1 )  =  1 

DO  16  J=1,NLAT 
16     MB(J+1)=2*J 

NLAT+2)=2*(NLAT+1) 

11=1 

18  J=2,NLONG 

MA(j}=2"(J-l) 
MAXNLONG+ 1 )  =  2*NLONG 
RETURN 
END 


14 


PROGRAM  -  POISSON'S 


APPENDIX  C 
EOUATION  with 
CONDITIONS 


DIRICHLET  BOUNDAR" 


Listing:   GAUSS 5  FORTRAN 


C 

C 

c 
c 

c 
c 
c 
c 
c 
c 
c 
c 


MAIN  PROGRAM 


STIFFNESS 
RESOLUTION 


MATRIX  USING  TENSOR  PRODUCT 


1000 


STIFFNESS  MATRIX 


TENSOR  PRODUCT 


503 
500 
1001 
C 

c 


501 

504 

1002 

1012 

1004 

1014 
C 


(24F3.0n 


NLAT  ='  13,/) 
AND  SA2  OF  STIFFNESS 


520 
521 
510 

C 
C 
C 

C 
C 


THIS  PROGRAM  IS  DESIGNED  TO  TEST  THE  SCHEME  WHICH 
RESOLVES  THE  STIFFNESS  MATRIX  INTO  A  SUM  OF  TWO  TENSOR 
PRODUCTS  IN  ORDER  TO  SOLVE  THE  SYSTEM  OF  EQUATIONS 
K  W  =  V.    THERE  ARE  DIRICHLET  BOUNDARYCONDITIONS  ON 
ALL  4  EDGES  OF  THE  REGION.   THE  PRESCRIBED  BOUNDARY 
VALUES  ARE  GIVEN  IN  THE  CORRESPONDING  LOCATIONS  IN  V. 
THE  SUBROUTINES  MAY  BE  INSERTED  IN  THE  PROGRAM  DEVISED 
BY  HINSMAN. 

IMPLICIT  REAL*8(A-H.O-Z) 
COMMON / CM1A/ NLAT , NLONG 
COMMON/ CM8 /A ( Z 1 ). B[Z 1] 

COMMON  AG(ZB],BG(ZC).GA1(ZK) ,SAl[ZK) .GBl(ZL) ,SB1(ZL) 
DIMENSION  V(ZP],Wl(z6) ,P(ZR) ,D(ZS) ,u(ZT) 
READ(5,  *) NLONG   NLAT 
LATX=NLAT+1 
WRITE (6, 1000) 
FORMAT?/ 7,  ' 
1RESOLUTION' ,/) 
WRITE? 6 , 100 1) NLONG , NLAT 
READ?  5. -O  A.  B 
WRITE (6, 5 00) A 
WRITEX6.503)B 
FORMAT?/,'    B. 
FORMAT (/,'    A:   '  {24F3.0 
FORMAT (      NLONG  =  ,13, ' 
CONSTRUCT  FACTORS,  GAl,  GB1,  SAl 
MATRIX 

CALL  AMTRX4 
WRITE (6, 501) AG 

V    AG:   ' 

504)BG 

BG '   ' 

1002)GA1 

3    GAl',/, 

l012)SAl 

,'   SAl',/, 

i004)GBl 

/.*   GB1',/ 

1614)SB1 
FORMAT (//,'   SB1'  / 
LOAD  BORDER  VECTOR  Wl 
READ(5.-)V 
CALL  BORDER (W1,V) 
WRITE(6.510)V 
Ll=4 "NLONG 
L2=L1+1 
L3=L1+4*LATX 

WRITE(6,520)(W1(L) ,L=1.L1) 
WRITEX6,521)(W1(L) ,L=L^.L3) 
FORMAT^/,    W1V,/  (3X,^F8.2)) 
FORMAT (3X5F8. 2) 

wi(L 

W1(L 

v>v,s 


FORMAT 
WRITE? 
FORMAT 
WRITE? 
FORMAT 
WRITE? 
FORMAT 
WRITE? 
FORMAT 
WRITE(6, 


■(12F4.1)) 
.(12F4.1)) 

(3X.6F7.3)) 
(3X,6F7.3)) 


(3X,6F7 
(3X,6F7 


3)) 
3)) 


CALL  MU 

WRITE (6 

WRITE 

WRITE 

CALL 

WRITE 

WRITE 

CALL 

WRITE 

WRITE 

WRITE 

READ( 


520 
6;521 
6,510 
ORDER 
6,520 
6,521 
ULT1H 
6,510 
6,520 
6,521 


5)) 


^L2,L3) 


,L=1. 

l,GBl) 


htt 


JL 


530 


)P 


W 


lJ!l=L2,L3) 


WRITE(6,530)P 

FORMAT(/,T       P:       ' , / , ( 6X , 3F12 . 4) ) 
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READC5  *)D 

WRITE?6,531)D 

531  FORMAT (/,'   D:   \3F12.4) 
C 

C   FORM  U  =  VINT*P 
C 

LONGM=NLONG-l 

LONGMM=NLONG-2 

NLATM=NLAT-1 

DO  30  L=1.NLATM 

JP1=(L-1)-NLATM 

KUl=CL-lj*(NL0NG-2)-l 

DO  29  K=2,LONGM 

TEMP=0. 

DO  28  J=1,NLATM 

JV=J-NLONG+K 

JP=JP1+J 

28  TEMP=TEMP+V(JV)*P(JP) 

29  U(KU1+K)=TEMP 

30  CONTINUE 
WRITE(6,532)U 

532  FORMAT??,1   U:   '  ,  /'  (6X,  4F12  .  4) ) 
CALL  FFFDB(U,D,GA1,SA1) 
WRITE(6,532)U 

C 

C   PUT  FINAL  RESULT  IN  V 

C 

DO  40  L=1,NLATM 

DO  3  9  K=l,LONGMM 

TEMP=0. 

DO  38  J=1,NLATM 

38  TEMP  =  TEMP  +  P((J-l)'vNLATM+L)'vU((J-l)»LONGMM+K) 

39  v7l*NLONG+K+1)=TEMP 

40  CONTINUE 
WRITE (6, 5 10  )V 
STOP 

END 
C 
C***********^**************^^ 

SUBROUTINE  BORDER (Wl.V) 

C 

C   THIS  SUBROUTINE  CLEARS  THE  BORDER  VECTOR  Wl  AND 

C   SUBSTITUTES  THE  BOUNDARY  VALUES  FROM  V. 

C 

IMPLICIT  REAL*8(A-H.O-Z) 

COMMON/ CM 1A/NLAT,NL0NG 

DIMENSION  W1(ZQ),V(ZP) 

LATX=NLAT+1 

NC=4*NLONG 

NR=4*LATX 

NB=NC+NR 

DO  4  J=1,NB 
4      W1(J)=0. 

DO  6  J=l,NLONG 
6      Wl(J)=vfj) 

DO  8  J=l,NLONG 

L=3"NLONG+J 

K=NLAT*NLONG+J 
8      W1(L}=V(K) 

DO  10  J=1,LATX 

L=NC+J 

K=(J-l)*NLONG+l 
10     W1(L)=V(K) 

DO  12  J=1,LATX 

L=NC+3*LATX+J 

K=J-NLONG 
12     W1?L)=V(K) 

RETURN 

END 

SUBROUTINE  AMTRX4 

C   THIS  SUBROUTINE  FORMS  THE  MATRICES  GA1,  GB1,  SA1,  AND  SB1 
C   THAT  ARE  FACTORS  IN  THE  TENSOR  PRODUCTS  USED  TO  FORM  THE 
C   COEFFICIENT  MATRIX  ("STIFFNESS"  MATRIX)  FOR  THE  POISSON 
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C   EQUATION.   ALL  OF  THESE  MATRICES  ARE  SYMMETRIC  AND 
C   TRIDIAGONAL. 


C 


IMPLICIT  REAL*8(A-H.O-Z) 
COMMON/ CM1A/NLAT , NLONG 
COMMON / CM^ / A  < Z 1 N<  B  ( Z 1 ) 


COMMON  AG(ZBl,B6lzC),GAl(ZK) , SAl(ZK) , GB1(; 
DIMENSION  BG(NLAT) , AG (NLONG) , GB1 (2*NLAT- I 


C 

C     lGAl(3*NLONG-3) 

C 

LATX=NLAT+1 

LONGM= NLONG -1 
C    FIND  BG  =  (ELEMENT  HEIGHT) /6 

NM= NLONG -1 

DO  2  J=1,NLAT 
2      BG(J)=B(i+NM*(J-l))/3. 
C    GENERATE  GB1  AND  6*SBi 

GBlfl)=2.*BG(l" 

SBlf  1)  =  1-  /BGU 

DO  4  J=2,NLAT 

K=2*(J-l) 

GB1(K)=2.*(BG(J-1)+BG(J)) 

GB1  K+1)=BG(J-1) 

SB1(K)=1./BG(J-1)+1./BG(J) 
4      SB1(K+1)=-1./BG(J-1) 

gb i  ( 2*nlat  )  =  2  .  *bg  ( nlat ) 
gbi(2"Nlat+i)=bg7nlat) 
sb1(2*nlat)=1. /bg(nlat) 

SBlf2*NLAT+l)=-l. /BG(NLAT) 

J2=2*NLAT+1 

DO  6  J=1,J2 
6      SBl(J)=SBl(J)/6. 
C    FIND  AG  =  (ELEMENT  WIDTH) /6. 

DO  10  J=l, LONGM 
10     AG(J)=A(j)/3. 
C    GENERATE  GA1  AND  6-SAl 

GA1(1)=2.*AG(1) 

SAlfl)=l./AG(l) 
DO  12  J= 2, LONGM 

K=2*(J-1) 

" (K)  =  2 


C  T!  "*  f  1 T   \ 


GAl(K)=2.*iAG(J-l)+AG(J)) 

GA1(K+1)=AG(J-1) 

SA1(K)=1./AG(J-1)+1./AG(J) 


GA1(K+1)=AG(J-1) 
SA1(K)=1./AG(J-1)+1. 
12     SAlfK+lt=-l./AG(J-l) 


GA1 ( 2 -L0NGM )  =  2 -AG (LONGM ) 
GA1 f 2 *LONGM+ 1 ) = AG (LONGM ) 
SA1(2'vL0NGM)  =  1./AG7L0NGM) 
SA1  f  2--LONGM+  1 )  =  -  SAI  ( 2 -LONGM) 
J2=2*NLONG-l 
DO  14  J=1,J2 
14     SAl(J)=SAi(J)/6. 
RETURN 
END 

C********  **  **   **       >.  *       *   ******  -**  ***• 

SUBROUTINE  MULT1 (Wl , V , A,B) 

C  SUBROUTINE  PREMULTIPLIES  Wl  MATRIX  BY  TRUNCATED  A  MATRIX 

C  (FIRST  AND  LAST  ROWS  OMITTED) ,  POSTMULTIPLIES  PRODUCT  BY 

C  TRUNCATED  B  MATRIX  7FIRST  AND  LAST  COLUMNS  OMITTED},  AND 

C  SUBTRACTS  INTERIOR  ELEMENTS  OF  Wl  FROM  CORRESPONDING 

C  ELEMENTS  OF  V. 

C 

IMPLICIT  REAL*8(A-H.O-Z) 

COMMON/CM  LA/ NLAT, NLONG 

DIMENSION  Wl(ZQ) ,V(ZP) ,A(1) ,B(1) 

LATX=NLAT+1 

LONGM=NLONG-l 

FCU=W1(1) 

RCU=W1(3-NL0NG+1) 

L1=4*NL0NG 

L2=L1+1 

L3=L1+4*LATX 

DO  2  J= 2, LONGM 

K=2*(J-l) 
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L=3*NL0NG+J 

FCC  =  A(K+ 1 ) *FCU  + A (K ) *W1 ( J ) + A (K+  3 ) *W1 ( J+ 1 ) 

RCC= A (K- 1 ) *RCU+A (K ) *W1 (L ) + A (K  + 3 ) *W1 (L+ 1 ) 

FCU=WI(JX 

RCU=W1(L 

W1(J)=FCC 
2      W1CL)=RCC 

LASTA=2*NL0NG-1 

DO  4  J=2,NLAT 

L=4*NL0NG+J 

LL=L+LATX 

K=L+3*LATX 

KK=K-LATX 

Wl(LL)=AC3)*Wl(L) 
4      Wl(KK)=AiLASTAj"Wl(K) 
C      WRITE 6,520)(Wl(L),L=l.Ll) 
C      WRITEI6.521)(W1(LJ ;l=L2,L3) 

520  FORMAT?//, '   INTERMEDIATE  RESULT,  Wl ' , / , ( 3X , 5F8 . 2 ) ) 

521  FORMAT (3X6F8. 2) 
LASTB=2*LATX-1 
NC=4*NL0NG 

Wl (NLONG+  2) =B ( 3 ) -Wl ( 2 ) +B (2 ) *W1 (NC+LATX+2 ) +B (5 ) 
1*WICNC+LATX+3J 

Wl ( 2 -NLONG- 1) =B ( 3 ) *W1 (NLONG -  1 ) +B ( 2 ) *W1 (NC  +  2 *LATX+  2 ) 
l+B(5)*Wl(NC+2*LATX+3) 

WlC2"NLONG+25=BfLASTB-2)*Wl(NC+2"LATX-2)+B(LASTB-3) 
1*WI7nC+2*LATX-  1  ]  +B(LASTB  ) -Wl  (  3 -'NLONG  +  2] 

Wlf3v-NLONG-l)=B(LASTB-2)"Wl(NC+3"LATX-2)+B(LASTB-3) 
1*WI (NC+3*LATX- 1 ) +B (LASTB ) *WI (NC- 1 ) 

J2=NLONG-2 

DO  6  J=3,J2 

W1(J  +  NL0NG>B(3)*W1(J) 
6      Wl  ( J  +  2 *NLONG ) =B  f  LASTB )  -Wl  ( J+  3*NLONG ) 

URL=Wl(NC+LATX+2) 

BRL=W1 [NC+2 "LATX+  2 ) 

J2=LATX-2 

DO  8  J=3,J2 

ND=2-(J-i) 

NCU=NC+LATX+J 

URC  =  B(ND+lV-vURL  +  B(ND)*Wl(NCU)+B(NDO)*Wl(NCU+l) 

BRC=B(ND+12"BRL+B(ND)"W1(NCU+LATX)+B(ND+3) 
1*W1(NCU+LATX+1) 

URL=W1(NCU) 

BRL=W1(NCU+LATX) 

W1(NCU)=URC 
8      W1(NCU+LATX}=BRC 

Wl ( NC  +  LATX+2 ) =  W1 (NLONG  +  2 ) 

Wl(NC+2"LATX-l)=Wl(2"NLONG+2) 

Wl C NC  + 2 -LATX+ 2 1 =W1 ( 2 "NLONG  -  1 1 

Wl (NC+  3*LATX- 1 ) =W1 ( 3 -NLONG- 1 ) 
C 

C   CORRECT  V 
C 

NLATM=NLAT-1 

DO  10  J=3,NLATM 

L=NC+LATX+J 

K=(J-l)"NLONG+2 

V(fc)=V(K)-Wl(L) 

L=LATX+L 

K=J*NLONG-l 

10    winm?1™ 

DO  12  J=2,J2 
L=NLONG+J 
V(L)=VTL)-W1(L) 
L=2-NL0NG+J 
K= (NLAT-l)*NLONG+J 
12     Y($)=V(K)-W1(L) 
RETURN 
END 

SUBROUTINE  FFFDB (X , E , GA, SA) 
C 

C   THIS  SUBROUTINE  SOLVES  A  SUCCESSION  OF  ONE-DIMENSIONSAL 
C   PROBLEMS.   THE  RELEVANT  COEFFICIENT  MATRIX  C  IS  FIRST 
C   FORMED,  THEN  FACTORED,  FOLLOWED  BY  FORWARD  REDUCTION, 
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C   DIVISION  BY  THE  DIAGONAL  ELEMENTS,  AND  BACK  SUBSTITUTION. 

C   THE  PROCESS  IS  CARRIED  OUT  NLATM  TIMES. 

C 

IMPLICIT  REAL*8(A-H.O-Z) 

COMMON/ CMIA / NLAT . NLONG 

DIMENSION  X(l) .E(l) ,GA(1) ,SA(1) ,C(ZU) 

NLATM = NLAT -t 

LONGM= NLONG -1 

LONGMM=NLONG-2 

DO  50  L=l, NLATM 
C 

C   FORM  COEFFICIENT* MATRIX  C 
C 

D1=E(L) 

C?1)=SA(2)+D1*GA(2) 

J2=2*NLONG-5 

DO  2  J=2,J2 
2      C(J)=SA(J+2)+Dl*GA(J+2) 
C 

C   FACTOR  C 
C 


TEMP=C(3)/C(1) 
C(2)=C(2i-TEMP* 

IF(C(2))7,7,3 


C?3)=TEMP 

J2=LONGMM-I 

DO  5  J=2,J2 

TEMP=C(2-J+1)/C(2*(J-1)) 


C 1 2 * J )  =  C ( 2 * J ) - TEMP"C ( 2*J+ 1 ) 

If7c(2*J))7,7 

Cf2*J+l)=TEMP 


7  WRITE(6,1000)J.C(2*J) 

1000   FORMAT  HZ  j.'  STOP  -  MATRIX  NOT  POSITIVE  DEFINITE'  //, 

1T  NONPOSITIVE  PIVOT  FOR  EQUATION   ',13,//,'  PIVOT  =  ', 
2D20.12) 
STOP 

C 

C   PERFORM  FORWARD  REDUCTION 

C 

8  J2=(L-l)"LONGMM 
DO  10  J=2.LONGMM 

10     X(J2+J)=x(j2+J)-X(J2+J-l)*C(2*(J-l)+l) 

C 

C   DIVIDE  BY  DIAGONAL  ELEMENTS 

C 

X(J2+1)=X(J2+1)/C(1) 

do  12  j=2:longmm 
12   x(j2+j)=x(j2+j)/c(2*(j-1)) 

C 

C   BACK- SUBSTITUTE 

C 

DO  14  J=2,LONGMM 

JB=J2+LONGM-J 
14     X( JB ) =x7 JB ) -X ( JB+ 1 ) *C ( 2 ~ (LONGM- J ) + 1 ) 
50     CONTINUE 

RETURN 

END 
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